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Abstract 

Image  segmentation  represents  an  essential  step 
in  the  early  stages  of  an  Automatic  Target 
Recognition  system.  We  propose  two  robust  ap¬ 
proaches  fundamentally  based  on  scale  informa¬ 
tion  inherent  to  a  given  imagery.  The  first  ap¬ 
proach  is  parametric  in  that  the  scale  evolution 
of  an  image  is  statistically  captured  by  a  model 
which  is  in  turn  utilized  to  classify  the  pixels 
in  the  image.  The  second,  on  the  other  hand, 
nonlinearly  evolves  the  image  along  some  spe¬ 
cific  characteristic  to  unravel  and  delineate  the 
various  comprising  entities  in  it. 

1  Introduction 

The  growing  interest  in  Automatic  Target  Re- 
congnition  (ATR)  is  primarily  due  to  its  impor¬ 
tance  in  many  applications  ranging  from  manu¬ 
facturing  to  remote  sensing  and  surface  surveil¬ 
lance.  Analysis  and  classification  of  various  en¬ 
tities  of  an  image  are  often  of  great  interest,  and 
its  partitioning  into  a  set  of  homogeneous  re¬ 
gions  or  objects  is  thus  of  importance.  The  mere 
size  of  some  imagery  constitutes  a  major  hur¬ 
dle.  A  case  in  point  is  Synthetic  Aperture  Radar 
(SAR)  imagery  for  which  terrain  coverage  rates 
are  very  high  (in  excess  of  1  km2 /s)  and  which 
with  daunting  computational  demands,  make 
algorithmic  efficiency  of  central  importance. 

A  SAR  image  reflects  a  coherent  integration  of 
scatterer  returns  (i.e.  reflectivity  characteris¬ 
tics)  within  a  resolution  cell.  The  number  of 
scatterers  which  coherently  sum  up  within  a 
cell  will  vary  with  the  resolution.  This  leads 
to  a  variation  in  the  underlying  statistics  of  the 
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image.  Natural  clutter  for  example,  tends  to 
consist  of  a  large  number  of  equivalued  scat¬ 
terers,  and  this  in  contrast  to  the  man-made 
one,  mostly  comprising  a  few  prominent  scat¬ 
terers.  It  is  precisely  this  type  of  statistical 
characteristic  that  we  are  interested  in  captur¬ 
ing  and  subsequently  using  as  the  basis  for  our 
classification  of  various  terrain  types  in  SAR. 
In  addressing  such  a  problem,  we  have  two  ba¬ 
sic  approaches,  the  first  aims  at  characterizing 
the  statistics  of  the  image  evolution  in  scale  and 
ultimately  using  them  in  the  pixel  classification; 
the  second  attempts  to  diffuse  “noncomplying” 
observations  via  a  nonlinear  evolution  to  make 
it  converge  to  specific  desired  domains  of  attrac¬ 
tion. 

The  goal  of  this  paper  is  to  address  a  funda¬ 
mental  problem  arizing  in  ATR,  namely  the  ro¬ 
bust  and  efficient  segmentation  of  imagery  into 
homogeneous  regions,  and  in  presence  of  per¬ 
haps  severe  noise.  To  proceed,  we  introduce  in 
the  next  section  a  multiscale  stochastic  model¬ 
ing  framework  which  affords  one  to  capture  the 
evolution  in  scale  of  a  given  image  process  and 
to  statistically  characterize  regions  of  interest. 
Two  algorithms  based  on  this  framework  will  be 
described  and  shown  to  lead  to  efficient  and  ac¬ 
curate  segmentation.  In  Section  4,  we  present 
our  second  method  based  on  a  variational  ap¬ 
proach,  and  which  consists  of  nonlinearly  evolv¬ 
ing  a  given  image  along  some  specific  geometric 
constraints  built  around  an  energy  functional. 
Finally,  in  Section  5,  we  provide  a  number  of  ex¬ 
amples  substantiating  the  proposed  algorithms 
using  real  data  imagery  and  accurate  estimates 
of  boudaries. 


2  Stochastic  Modeling 

2.1  Multiscale  Stochastic  Models 

In  this  subsection,  we  describe  a  general  multi¬ 
scale  modeling  framework  e.g.  [2]  and  its  adap¬ 
tation  to  classification/identification  problems. 
Under  this  framework,  a  multiscale  process  is 
mapped  onto  nodes  of  a  qth  order  tree,  where 
q  depends  upon  how  the  process  progresses  in 
scale. 


As  illustrated  in  Fig.  1,  a  qth  order  tree  is  a 
connected  graph  in  which  each  node ,  starting  at 
some  root  node,  branches  off  to  q  child  nodes. 
As  described  above,  the  appropriate  represen¬ 
tation  for  a  multiscale  SAR  image  sequence  is 
q  =  4,  a  quadtree.  Each  level  of  the  tree  (i.e., 
distance  in  nodes  from  the  root  node)  can  be 
viewed  as  a  distinct  scale  representation  of  a 
random  process,  with  the  resolutions  proceeding 
from  coarse  to  fine  as  the  tree  is  traversed  from 
top  to  bottom  (root  node  to  terminal  nodes).  A 
coarse-scale  shift  operator,  7  is  defined  to  ref¬ 
erence  the  parent  of  node  s,  just  as  the  shift 
operator  z  allows  referencing  of  previous  states 
in  discrete  time-series.  The  state  elements  at 
these  nodes  may  be  modeled  by  the  coarse-to- 
fine  recursion 

x(s)  =  A(s)x(s7)  +  B(s)w(s).  (1) 


at  any  node  s,  the  processes  defined  on  each  of 
the  distinct  subtrees  extending  away  from  node 
s  are  mutually  independent. 

For  pixel  classification  purposes,  a  multiscale 
model  can  be  constructed  for  each  specific  ho¬ 
mogeneous  class.  To  specify  each  model,  it  is 
necessary  to  determine  the  appropriate  coeffi¬ 
cients  in  the  matrices,  A(s)  and  B(s),  and  the 
statistical  properties  of  the  driving  noise,  w(s). 
Once  the  models  have  been  specified,  a  likeli¬ 
hood  ratio  test  can  be  derived  to  segment  the 
imagery  into  the  clutter  classes. 

For  a  binary  classification  problem  (i.e.  requir¬ 
ing  a  binary  hypthesis  test)  each  pixel  in  the  im¬ 
age  corresponds  to  one  of  two  hypotheses:  the 
pixel  is  part  of  some  texture  ( Hg )  or  another 
(Hf).  By  exploiting  the  Markov  property  asso¬ 
ciated  with  the  multiscale  models  for  imagery, 
the  log-likelihood  ratio  test  for  classifying  each 
pixel  can  be  written  as, 


1  =  E,  loS  |Px(*)|/39  >  (x(fi)  I  0g) 
£slog  [px(s)i/3/(X(s)  I  (3f) 


(2) 


where  /3(.}  =  (X(sj ),#(.)),  and  px(s)\0g  and 
Px(s)\/3j  are  the  conditional  distributions  for  the 
two  hypothesized  models.  In  the  next  subsec¬ 
tion,  we  will  show  that  this  likelihood  test  can 
be  efficiently  computed  in  terms  of  the  distribu¬ 
tions  for  w(s)  under  the  two  hypotheses. 


2.2  Scale-Autoregressive  SAR  Model 


In  this  paper,  we  focus  on  a  specific  class  of 
multiscale  models,  namely  scale-autoregressive 
models  [2]  of  the  form 


In  this  recursion,  A(s)  and  B(s)  are  matrices  of 
appropriate  dimension  and  the  term  w(s)  repre¬ 
sents  white  driving  noise.  The  matrix  A(s)  cap¬ 
tures  the  deterministic  progression  from  node 
sq  to  node  s,  i.e.,  the  part  of  x(s)  predictable 
from  x(s7),  while  the  term  B(s)w(s)  represents 
the  unpredictable  component  added  in  the  pro¬ 
gression.  An  attractive  feature  of  this  frame¬ 
work  is  the  efficiency  it  provides  for  signal  pro¬ 
cessing  algorithms.  This  stems  from  the  Markov 
property  of  the  multiscale  model  class,  which 
states  that,  conditioned  on  the  value  of  the  state 


I(s)  =  ai(s)I(s7)  +  a2(s)I(s^f2)  +  . . .  + 

aR(s)I{s7yR )  +  w(s),  a{(s)  G  JR  (3) 

where  w(s)  is  white  driving  noise  and  “s”  is 
a  three-tuple  vector  denoting  scale  (or  resolu¬ 
tion  level),  and  spatial  coordinates  ( q,n ).  For 
homogeneous  regions  of  texture,  the  predic¬ 
tion  coefficients  (the  al(s)  in  (3))  are  constant 
with  respect  to  image  location  for  any  giver- 
scale.  That  is,  the  coefficients,  ai(s), . . .  ,  an(s), 
depend  only  on  the  scale  of  node  s  (de¬ 
noted  by  m(s)),  and  thus  will  be  denoted  by 


ai,ro(s)>  •  •  •  1  aR,m(s)-  Furthermore,  the  probabil¬ 
ity  distribution  for  w(s)  depends  only  on  m(s). 
Thus,  specifying  both  the  scale-regression  coef¬ 
ficients  and  the  probability  distribution  for  w(s ) 
at  each  scale  completely  specify  the  model. 

Following  the  procedure  of  state  augmentation 
used  in  converting  autoregressive  time  series 
models  to  state  space  models,  we  associate  to 
each  node  s  a  R-dimensional  vector  of  pixel  val¬ 
ues,  where  R  is  the  order  of  the  regression  in 
(3).  The  components  of  this  vector  correspond 
to  the  SAR  image  pixel  associated  with  node  s 
and  its  first  R  —  1  ancestors.  Specifically,  we 
define 

x(s)  =  [  I(s)  I(s7)  •••  I{sjR~l)  ]T.  (4) 

Thus,  for  a  model  of  the  form  (3)  or  equivalently 
(1),  i  in  (2)  can  be  calculated  using 

Px(s)|x(s7) (X(s)  |  X(sy))  —  Poi(s)(^'F(s))(5) 

with  W(s)  being  the  vector  of  residuals  at  scale 
“s”.  The  identification  of  the  model  for  each 
clutter  class  can  thus  be  obtained  for  each  scale 
m  by  a  standard  least-squares  minimization, 


ar 


where 


E  u«- 

{s  |  m(s)=m} 

-  ...  -  aRimI(sjR)]2^Q) 


arg  mm 
ameMR 


a  i 


am  —  [al,m  a2,m  ■■■  aR,m\  i 

with  R  being  the  regression  model  order.  For 
most  of  the  results  presented  in  Section  5  a  third 
order  regression  {R  =  3)  for  both  the  grass  and 
the  forest  models  was  chosen.  To  obtain  a  sta¬ 
tistical  characterization  of  the  prediction  error 
residuals  (the  w(s)  in  (3))  of  the  model  at  scale 
m,  we  evaluate  the  residuals  in  predicting  scale 
m  of  a  homogeneous  test  region.  In  particu¬ 
lar,  we  use  the  am,;i  found  in  (6)  to  evaluate  all 
{w(s)|m(s)  =m}  as  indicated  by  Eq.  (5)  with 
y,  specifying  “grass”  or  “forest” . 


3  Model-Based  Classification 


3.1  Residual-Based  Classification 


While  we  could  conceivably  postulate  a  spatial 
random  field  model  for  each  windowed  clutter 


category  to  classify  its  center  pixel,  we  use  to 
advantage  the  efficiency  of  multiscale  likelihood 
calculation  to  base  the  classification  of  each  in¬ 
dividual  pixel  “s”  on  a  surrounding  (2 K  +  1)  x 
(2 K  +  1)  window  W(s),  where  the  parameter 
AT  is  a  nonnegative  and  judiciously  selected  in¬ 
teger.  While  a  larger  window  provides  a  more 
accurate  classification  of  homogeneous  regions, 
it  also  increases  the  likelihood  that  the  window 
contain  a  clutter  boundary.  Thus,  keeping  the 
window  size  as  small  as  possible  is  also  desir¬ 
able.  As  shown  in  [1]  in  the  next  section,  one 
can  determine  the  tradeoff  between  classifica¬ 
tion  accuracy  and  window  size  by  examining  the 
empirical  distribution  of  l  over  windows  of  var¬ 
ious  size  for  homogeneous  regions  of  grass  and 
forest. 


3.1.1  Statisical  Hypotheses  Test 


Whenever  a  clutter  boundary  is  present  within  a 
test  window  the  validity  of  the  center  pixel  clas¬ 
sification  is  questionable.  This  effect  results  in  a 
classification  bias  near  boundaries.  To  address 
this  problem,  we  devise  a  method  to  detect  the 
proximity  of  grass-forest  boundaries  and  subse¬ 
quently  utilize  a  procedure  to  refine  the  classifi¬ 
cation.  Terrain  boundary  proximity  is  detected 
via  a  simple  modification  of  the  decision  made 
based  on  the  test  statistic  t.  Specifically,  rather 
than  comparing  l  to  a  single  threshold  to  decide 
on  a  grass-or-forest  classification,  we  compare  £ 
to  two  thresholds  a  and  b  as  shown  in  Fig.  2, 
and  where  we  include  a  defer  decision. 


Fig.  2:  Initial  pixel  classification. 

This  is  tantamount  to  refining  the  decision  pro¬ 
cedure  by  considering  smaller  windows  within 
the  original  window  until  a  majority  rule  lifts 
the  ambiguity, 


Hg  :  l  >  a  Classify  as  Grass. 

a  >  t  >  b  Defer  decision 

(possible  boundary  presence), 


Hf  :  £  <  b  Classify  as  Forest. 


3.2  Parameter-Based  Classification 

An  alternative  approach  to  classifying  a  clut¬ 
ter  type  is  to  evaluate  an  l2  distance  between 
the  computed  model  parameters  for  a  window 
region  W(s)  and  those  of  a  template  homoge¬ 
neous  region  [3].  For  each  pixel  “s”  (or  equiv¬ 
alently  tree  node)  a  model  for  a  corresponding 
surrounding  window  W(s)  is  computed  at  sev¬ 
eral  scales  giving  rise  to  what  we  refer  to  as  an 
evolution  vector  statistic, 

y(®)  [^*m(s)>  C*m(s)  —  l  i  '  '  ’  J  qu];  (7) 

where  “s”  will  sweep  all  of  the  three-tuple  vec¬ 
tors  [m(s),g,  n]  at  the  m(s)th  resolution  with 
m{s)  —  I,--  -  ,L  —  1  and  L  is  the  cardinality 
of  the  considered  set  of  images.  We  should 
note  at  this  point  that  in  this  approach,  the 
DC  component  in  the  is  included.  The  or¬ 
der  of  the  regression  associated  with  modeling 
W([m(s),  q,  n])  from  its  ancestors  will  vary  with 
the  level  m(s)  as  defined  by  the  function 

Om(s)  =  ma x(i?,m(s)) 

The  regression  vector  a(m(s))  provides  a  sta¬ 
tistically  optimal  description  of  the  linear  de¬ 
pendency  of  W(s)  on  {W(s7J)}J0=T’,  and  Y(s) 
thus  being  a  measure  of  the  scale  evolution  be¬ 
havior  of  the  windowed  region. 

Here  again  the  size  of  the  window  used  to  com¬ 
pute  the  evolution  vector,  is  the  result  of  a 
tradeoff  between  modeling  consistency  and  lo¬ 
cal  accuracy. 

3.2.1  Statistical  Classification 

A  characterization  of  the  evolution  vector  y(s) 
is  necessary  to  carry  out  a  statistically  mean¬ 
ingful  classification  of  various  types  of  terrain. 
Specifically,  a  BHT  is  applied  to  the  evolution 


y(s)  in  order  to  classify  node  s  as  a  member 
of  either  a  region  of  grass  or  of  forest.  These 
hypotheses  which  are  respectively  designated 
as  hypotheses  Hg  and  Hf.  The  classification 
of  pixel  m(s)  will  depend  only  on  y(s)  and 
the  predetermined  likelihoods  p(y(s)\Hg)  and 
p(y(s)\Hf). 

To  carry  out  a  statistically  significant  hypoth¬ 
esis  test,  we  need  to  specify  the  conditional 
probability  density  for  y  under  each  hypothe¬ 
sis.  To  do  this,  we  extensively  examined  the 
distribution  of  the  evolution  vectors  obtained 
from  a  large  homogeneous  region  of  the  cor¬ 
responding  terrain.  For  both  grass  and  forest 
terrain,  it  turns  out  that  each  component  in 
y  approximately  has  a  Gaussian  distribution. 
We  consequently  make  the  approximation  that 
p(y\Hg)  and  p  (y\Hf)  are  N- variate  Gaussian 
densities.  They  are  then  completely  specified  by 
their  mean  vectors  m5  and  m j  and  their  covari¬ 
ance  matrices  Kg  and  K f  all  of  which  are  calcu¬ 
lated  from  the  training  data  for  each  hypothesis. 
A  maximum  likelihood  (ML)  detector  is  then 
used  to  classify  each  pixel  (using  the  ML  detec¬ 
tor  assumes  equal  apriori  probabilities  for  each 
hypothesis  and  a  cost  function  that  penalizes  all 
misclassifications  equally) .  In  implementing  the 
ML  detector,  a  threshold  p  is  calculated  from 
the  likelihoods  and  used  in  the  classification  of 
each  pixel  through  its  comparison  to  a  sufficient 
statistic  derived  from  the  evolution  vector.  Be¬ 
cause  p(y\Hg)  and  p(y\Hf)  are  assumed  to  be 
Gaussian,  it  is  straightforward  to  compute  the 
threshold  p  and  the  sufficient  statistic  i'  (y)  for 
each  evolution  vector  as 

r/==1°S}^|  ’ 

and 

?  (y)  =  (y  -  m/)T  KJ1  (y  -  mf)  - 

(y  —  ms)T  K~l  (y  —  ms) .  (9) 

The  classification  of  I(s),  denoted  as  C(s),  is 
then  given  by 


P 


Hg  if  P  >  i  (y [,,„]) 

Hf  iip<i  (y [9,n]) . 


(10) 


The  construction  of  the  evolution  y  and  sub¬ 
sequent  application  of  a  BHT  for  each  s  € 


SIDE 


{s  |  m(s)  =  m}  thus  provide  a  segmentation 
of  Im.  Instead  of  independently  generating  a 
segmentation  for  each  image  resolution  for  all 
s  6  {s  j  m(s )  <  L},  C(s)  can  be  obtained 
by  comparing  rj  to  the  average  of  the  sufficient 
statistics  of  nodes  in  I/,  which  have  s  as  an  an¬ 
cestor,  i.e. 


v  >  5^(y(s)) 

S 

v  <  J^(y(s)). 

s 

(ii) 


Although  this  does  not  yield  a  sufficient  statis¬ 
tic  for  I(s),  doing  so  is  computationally  more 
efficient  than  calculating  a  sufficient  statistic  as 
in  Eq.  (9)  for  each  node  at  every  level  in  the 
quadtree.  Note  that  the  segmentation  technique 
described  here,  could  easily  be  generalized  to  a 
larger  number  of  terrain  types. 


4  Stabilized  Inverse  Diffusion  Equa¬ 
tions  (SIDEs) 


The  previous  two  techniques  rely  on  modeling  a 
linear  evolution  of  the  observed  imagery  with  an 
ultimate  goal  of  pixel  classification.  In  this  sec¬ 
tion,  we  instead  carry  out  a  nonlinear  evolution 
which  is  driven  by  prescribed  geometric  features 
underlying  the  process/imagery,  and  study  its 
progression. 

Towards  that  end,  we  introduce  a  discontinuous 
force  function,  resulting  in  a  system  of  equations 
that  has  discontinuous  right-hand  side  (RHS). 
As  shown  below,  the  objective  is  to  drive  an 
evolution  trajectory  onto  a  lower-dimensional 
surface  which  clearly  has  value  in  image  anal¬ 
ysis,  and  in  particular  in  image  segmentation. 
Segmenting  a  signal  or  image,  represented  as  a 
high-dimensional  vector  I,  consists  of  evolving 
it  so  that  it  is  driven  onto  a  comparatively  low¬ 
dimensional  subspace,  which  corresponds  to  a 
segmentation  of  the  signal  or  image  domain  into 
a  small  number  of  regions. 


Pcrona-MaJik 
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Fig  3:  Force  function  for  a  SIDE. 


The  type  of  force  function  of  interest  to  us  here 
is  illustrated  in  Fig.  3 (right).  More  precisely, 
we  wish  to  consider  force  functions  F(v)  which, 
in  addition  to  drving  the  following  evolutions, 


i  (s)  =  F(I)(s),  (12) 

1(0)  =  I0, 

where  “s”  is  now  a  continuous  scale,  satisfy  the 
following  properties: 

F'(v)  <  0  foru  ^  0,  and  i?(0+)  >  0, 

F( vi)  =  F(v 2)  <£>  v\  =  V2-  (13) 

Contrasting  this  form  of  a  force  function  to  the 
Perona-Malik  function  [4]  in  Fig.  3  (left),  we 
see  that  in  a  sense,  one  can  view  the  discon¬ 
tinuous  force  function  as  a  limiting  form  of  the 
continuous  force  function  in  Fig.  3  (left).  We 
therefore  need  a  special  definition  of  how  the 
trajectory  of  our  evolution  proceeds  at  the  dis¬ 
continuity  points  F(0+)  ^  F(0~)).  For  this  def¬ 
inition  to  be  useful,  the  resulting  evolution  must 
satisfy  well-posedness  properties:  the  existence 
and  uniqueness  of  solutions,  as  well  as  stabil¬ 
ity  of  solutions  with  respect  to  the  initial  data. 
Assuming  the  resulting  evolutions  to  be  well- 
posed,  we  demonstrate  that  they  have  the  qual¬ 
itative  properties  we  desire,  namely  that  they 
both  are  stable  and  also  act  as  inverse  diffu¬ 
sions  and  hence  enhance  edges.  We  address  the 
issue  of  well-posedness  and  other  properties  in 
[5]. 

Considering  the  evolution  (12)  with  F(v)  as  in 
Fig.  3  (right)  in  a  SIDE,  notice  that  the  RHS  of 
(12)  has  a  discontinuity  at  a  point  I  if  and  only 
if  F  —  I{+ 1  for  some  i  between  1  and  N  —  1. 
It  is  when  a  trajectory  reaches  such  a  point  I(.) 
that  we  need  the  following  definition: 

if  =  ii+i  =  h(r(i.+2  -  ii+0  -  F(ii  - 1.-1)). 

(14) 


In  other  words,  the  two  observations  are  simply 
merged  into  a  single  one,  resulting  in  Eq.  14  for 
n  =  i  and  n  =  i  +  1  (the  differential  equations 
for  n  7^  i  +  1  do  not  change.). 

Similarly,  if  q  consecutive  observations  become 
equal,  they  are  merged  into  one  which  is 
weighted  by  1/q  [5].  The  evolution  can  then 
naturally  be  thought  of  as  a  sequence  of  stages: 
during  each  stage,  the  right-hand  side  of  (12) 
is  continuous.  Once  the  solution  hits  a  discon¬ 
tinuity  surface  of  the  right-hand  side,  the  state 
reduction  and  re-assignment  of  qni's,  described 
above,  takes  place.  The  solution  then  proceeds 
according  to  the  modified  equation  until  the 
next  discontinuity  surface,  etc. 

Notice  that  such  an  evolution  automatically 
produces  a  multiscale  segmentation  of  the  origi¬ 
nal  signal  if  we  view  each  compound  observation 
as  a  region  of  the  signal.  The  algorithm  may  be 
be  summarized  as  follows: 

1.  Start  with  the  trivial  initial  segmentation: 
each  sample  is  a  distinct  region. 

2.  Evolve  (12)  until  the  values  in  two  or  more 
neighboring  regions  become  equal. 

3.  Merge  the  neighboring  regions  whose  val¬ 
ues  are  equal. 

4.  Go  to  step  2. 

The  same  algorithm  can  be  used  for  2-D 
images,  which  is  immediate  upon  re-writing 
Eq.  (12)  and  an  example  is  provided  next: 

5  Experiments 

A  number  of  segmentation  experiments  have 
been  carried  out  on  SAR  imagery  using  the 
three  techniques  described  above.  Fig.  4  shows 
a  segmentation  as  a  result  of  a  systematic  pixel 
classification  as  described  in  Technique  1.  A 
brute  force  as  well  as  a  refined  segmentations 
are  shown,  and  the  handrawn  of  the  eyeballed 
boundary  is  shown  in  white.  In  Fig.  5,  the 
model-based  segmentation  is  illastrated  and  in 
Fig.  6  the  nonlinear  evolution-based  segmenta¬ 
tion  is  shown  which  clearly  is  perhaps  the  most 


promising,  albeit  at  a  slightly  higher  computa¬ 
tional  cost. 


Fig.  4:  Residual  and  Model-Based  Segmentation 
(b  and  c). 


Fig.  6:  Segmentation  by  Nonlinear  Diffusion. 
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